# Import libraries
library(plotly)
library(heatmaply)
library(chromVAR)
library(motifmatchr)
library(BSgenome.Hsapiens.UCSC.hg19)
library(SummarizedExperiment)
library(BuenColors)
library(chromVARxx)
library(chromVARmotifs)
library(data.table)
library(BiocParallel)
register(MulticoreParam(4)) # adjust according to your machine
peakfile <- "../../data/fromMA/peak-ids.bed"
peaks <- getPeaks(peakfile)
counts_ma_in <- fread("../../data/fromMA/combined-counts-clean-subject-corrected-quantiles-ma.tsv")
metadata_ma <- fread("../../data/fromMA/metadata-clean-ma.tsv")
counts_mat <- round(data.matrix(data.frame(counts_ma_in)[,-60]),0)
# Make Summarized Experiment
counts <- SummarizedExperiment(
rowRanges = peaks,
colData = metadata_ma,
assays = list(counts = counts_mat)
)
counts <- addGCBias(counts, genome = BSgenome.Hsapiens.UCSC.hg19)
data("human_pwms_v2") # also human_pwms_v1 with more motifs (typically redunant)
motif_ix <- matchMotifs(human_pwms_v2, counts, genome = BSgenome.Hsapiens.UCSC.hg19)
dev <- computeDeviations(object = counts, annotations = motif_ix)
variabilityAll <- computeVariability(dev)
plotVariability(variabilityAll, use_plotly = TRUE)
Figure 1a; variability across motifs
mostvariable <- tail(sort(variabilityAll$variability, index.return = TRUE)$ix,30)
m <- assays(dev)[["z"]][mostvariable,]
rownames(m) <- rowData(dev)$name[mostvariable]
heatmaply(m, colors = jdb_palette("solar_extra"))
Figure 2a; motifs x samples
#+ cache = FALSE, message=FALSE, warning=FALSE, fig.height = 10, fig.width = 10, echo = TRUE, fig.align=‘center’, fig.cap = “Figure 3a; sample tSNE in motif space”
plotDeviationsTsne(dev, deviationsTsne(dev, perplexity = 10, threshold = 2.5), sample_column = "cell.type",
shiny = FALSE)[[1]] %>% plotly::ggplotly()